###Replication Code for field experimental results in "Competence versus Priorities: Negative Electoral Responsesto Education Quality in Brazil"
###Authors: Taylor Boas, F. Daniel Hidalgo, Guillermo Toral
###Date: June 23, 2020
#R Version, Package Versions, and Session Information Located at the Bottom of File

#Set Working Directory to directory containing this file and "data" subdirectory. Directory must also have  "tables" and "figures" subdirectories.

#Following packages need to be installed: "tidyverse"; "cowplot"; "huxtable"; "estimatr", "future", "texreg", " "glmnet", "huxtable", "recipes", "kableExtra", "parsnip", "yardstick", "rsample", "furrr"


##Setup
set.seed(1202019)

library("tidyverse"); library("cowplot"); library("huxtable")
library("estimatr"); library("texreg"); library("kableExtra")

source("lasso_covars.R")

demean <- function(x){
  x - mean(x, na.rm = TRUE)
}


# Import Data and Create Estimation Datasets ------------------------------

panel <- read_csv("./data/panel_cleaned.csv")

#Create ipw weights
panel <- group_by(panel, codesetor) %>%
  mutate(prob_treat = mean(ana),
         ipw = ifelse(ana == 1, 1 / prob_treat, 1 / (1-prob_treat))) %>%
  ungroup()

panel_noattrit <- panel

covars <- c("age", "politics_interest", "turnout_2012",
            "vote_2012", "turnout_2014",
            "vote_2014", "partisan", "muni_biggest_prob", "politician_helped",
            "govt_eval_baseline", "acc_eval_baseline",
            "uncertain_acc_baseline", "edu_eval_baseline",
            "uncertain_edu_baseline", "tce_knowledge", "ana_knowledge",
            "child_school", "confid_fedgov",
            "confid_justice", "confid_tce", "confid_muni",
            "acc_responsible", "edu_responsible",
            "prob_vote_count", "prob_vote_monitoring",
            "prob_vote_count", "acc_rejected_prior", "tce_prior_cert",
            "ana_prior", "edu_prior_uncert",
            "years_edu", "race", "religion", "income", "relative_wellbeing",
            "female", "edu_rank1", "acc_rank1")

num_perms <- 2000

panel <- filter(panel, attrited != 1)
##Impute missing data for pretreatment covariates
covar_data <- mlr::impute(panel[,covars],
                         classes = list(integer = mlr::imputeMean(), numeric = mlr::imputeMean(),
                                        factor = mlr::imputeMode(), character = mlr::imputeMode()))$data[, covars]

panel <- bind_cols(panel[, ! names(panel) %in% covars], covar_data)
panel$codesetor <- as.factor(panel$codesetor)

##Create gap between reality and prior
panel$ana_prior_gap <- panel$ana_rank - panel$ana_prior

#Create tercile bins for performance metric
panel <- mutate(panel, ana_tercile = fct_recode(as.character(ntile(ana_rank, 3)),
                                              "3rd Tercile" = "3",
                                              "2nd Tercile" = "2",
                                              "1st Tercile" = "1"))
##Create tercile bins for gaps between performance and prior
panel <- mutate(panel, priorgap_tercile = fct_recode(as.character(ntile(ana_prior_gap, 3)),
                                                     "3rd Tercile" = "3",
                                                     "2nd Tercile" = "2",
                                                     "1st Tercile" = "1"))
##Create subsets of the full panel dataset
ana_contr_all <- filter(panel,  accounts != 1 &
                               attrited != 1)
ana_contr_good <- filter(panel, good_ana == 1 &
                               accounts != 1 & attrited != 1)
ana_contr_bad <- filter(panel, good_ana == 0 &
                              accounts != 1 & attrited != 1)

##Recipe package function to create estimation dataset
make_estimdata <- function(data, dv, covars){
  data <- data[, names(data) %in% c(covars, dv, "ana", "ipw", "ana_tercile", "ana_rank")]
  rec <- recipe(as.formula(paste(dv, " ~ ", paste(c(covars), collapse = "+"))),
                      data = data) %>%
    step_dummy(all_nominal()) %>%
    step_meanimpute(all_predictors()) %>%
    step_zv(all_predictors()) %>%
    prep()
  estim_data <- juice(rec)
  estim_data$ana <- data$ana
  estim_data$ana_rank <- data$ana_rank
  estim_data$ana_tercile <- data$ana_tercile
  estim_data$ipw <- data$ipw
  estim_data
}


##Use this function to estimate effects within tercile bins
est_func <- function(df){
  covars <- lasso_covars(data = df, dv = "vote",
                             covars = c(covars, "codesetor"), treat = "ana")
  results <- lm_lin(vote ~ ana, covariates = as.formula(paste("~",
                                                          paste(covars, collapse = " + "))),
                data = df, se_type = "HC1", weights = ipw) %>%
    tidy() %>%
    filter(term == "ana")
  results
}



# Distribution of ANA Performance Change ----------------------------------

load("./data/pe_munis.RData")

ggplot(pe_munis, aes(x = read.math.change)) +
  geom_histogram() +
  theme_minimal() +
  xlab("Change in Average ANA Test Score")
ggsave("./figures/ana_hist.pdf", width = 6, height = 4)

# Effect by ANA RANK ------------------------------------------------------

edu_estimdata <- make_estimdata(ana_contr_all, dv = "vote", covars = c(covars, "codesetor"))
edu_covars <- lasso_covars(data = edu_estimdata, dv = "vote",
                           covars = c(covars, "codesetor"), treat = "ana")
edu <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank +",
                                      paste(edu_covars, collapse = " + "))),
       data = edu_estimdata, weights = ipw)

edu_nocovars <- lm_robust(vote ~ ana * demean(ana_rank),
                          data = ana_contr_all, weights = ipw,
                          fixed_effects = codesetor)

##Divide Sample into Tercile by Ana Rank and estimate effect in each bin
bytercile_rank <- group_by(edu_estimdata, ana_tercile) %>%
  nest()
bytercile_rank$median <- map_dbl(bytercile_rank$data, ~ median(.$ana_rank))
bytercile_rank <- bytercile_rank %>%
  mutate(results = map(data, est_func))
bytercile_rank <- unnest(bytercile_rank, results)


##Calculate Marginal Effects Using Continuous Rank
get_ame <- function(data, mod){
  ame <- coef(mod)["ana"] + (data$ana_rank - mod$scaled_center["ana_rank"]) * coef(mod)["ana:ana_rank_c"]
  ame_var <- vcov(mod)["ana", "ana"] +
    (data$ana_rank - mod$scaled_center["ana_rank"])^2 * vcov(mod)["ana:ana_rank_c","ana:ana_rank_c"] +
    2 * (data$ana_rank - mod$scaled_center["ana_rank"]) * vcov(mod)["ana:ana_rank_c", "ana"]
  ame_plotdata <- tibble(mod = data$ana_rank, ame = unlist(ame), ame_se = sqrt(unlist(ame_var)))
  ame_plotdata
}

edu_ame <- get_ame(edu_estimdata, edu)

edu_ameplot <- ggplot(edu_ame, aes(x = mod, y = ame)) +
  geom_line(size = 1) +
  geom_hline(yintercept = 0, lty = 2, size = 1) +
  geom_ribbon(aes(ymin = ame - 1.96 * ame_se,
                  ymax = ame + 1.96 * ame_se), alpha = .1) +
  geom_pointrange(data = bytercile_rank, aes(x = median, y = estimate,
                                             ymax = estimate + 1.96 * std.error,
                                             ymin = estimate - 1.96 * std.error),
                  size = .5) +
  theme_minimal() +
  ylab("Average Marginal Effect of Information") +
  xlab("Standardized Test (ANA) Performance Rank")
edu_ameplot
xhist <- axis_canvas(edu_ameplot, axis = "x") +
  geom_histogram(data = edu_ame, aes(x = mod), binwidth = 5)
edu_ameplot <- insert_xaxis_grob(edu_ameplot, xhist, position = "bottom")

pdf("./figures/ame_rank.pdf", width = 6, height = 4)
ggdraw(edu_ameplot)
dev.off()

#EPS Version
cairo_ps(filename = "./figures/ame_rank.eps",
         width = 6, height = 4, pointsize = 12,
         fallback_resolution = 300)
ggdraw(edu_ameplot)
dev.off()



# Effects by Child in School ----------------------------------------------


edu_child_estimdata <- make_estimdata(filter(ana_contr_all, child_school ==1),
                                      dv = "vote", covars = c(covars, "codesetor"))
edu_child_covars <- lasso_covars(data = edu_child_estimdata,
                                 dv = "vote",  covars = c(covars, "codesetor"), treat = "ana")
edu_child <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank +",
                                                        paste(edu_child_covars, collapse = " + "))),
              data = edu_child_estimdata, weights = ipw)

edu_child_nocovars <- lm_robust(vote ~ ana * demean(ana_rank),
                             data = filter(ana_contr_all, child_school == 1),
                             fixed_effects = codesetor, weights = ipw)

edu_nochild_estimdata <- make_estimdata(filter(ana_contr_all, child_school == 0),
                                      dv = "vote", covars = c(covars, "codesetor"))
edu_nochild_covars <- lasso_covars(data = edu_nochild_estimdata,
                                 dv = "vote",  covars = c(covars, "codesetor"), treat = "ana")
edu_nochild <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank +",
                                                              paste(edu_nochild_covars, collapse = " + "))),
                    data = edu_nochild_estimdata, weights = ipw)


edu_nochild_nocovars <- lm_robust(vote ~ ana * demean(ana_rank),
                                data = filter(ana_contr_all, child_school == 0),
                                fixed_effects = codesetor, weights = ipw)

bytercile_child <- group_by(edu_child_estimdata, ana_tercile) %>%
  nest()
bytercile_child$median <- map_dbl(bytercile_child$data, ~ median(.$ana_rank))
bytercile_child <- bytercile_child %>%
  mutate(results = map(data, est_func))
bytercile_child <- unnest(bytercile_child, results)


bytercile_nochild <- group_by(edu_nochild_estimdata, ana_tercile) %>%
  nest()
bytercile_nochild$median <- map_dbl(bytercile_nochild$data, ~ median(.$ana_rank))
bytercile_nochild <- bytercile_nochild %>%
  mutate(results = map(data, est_func))
bytercile_nochild <- unnest(bytercile_nochild, results)


##Estimate Linear Interaction for Respondents with Children
child_ame <- get_ame(edu_child_estimdata, edu_child)
##Estimate Linear Interaction for Respondents without Children
nochild_ame <- get_ame(edu_nochild_estimdata, edu_nochild)

child_ameplot <- ggplot(child_ame, aes(x = mod, y = ame)) +
  geom_line(size = 2) +
  geom_hline(yintercept = 0, lty = 2, size = 1) +
  geom_ribbon(aes(ymin = ame - 1.96 * ame_se,
                  ymax = ame + 1.96 * ame_se), alpha = .1) +
  geom_pointrange(data = bytercile_child, aes(x = median, y = estimate,
                                        ymax = estimate + 1.96 * std.error,
                                        ymin = estimate - 1.96 * std.error),
                  size = .5) +
  theme_minimal() +
  ylab("Average Marginal Effect of Information") +
  xlab("Standardized Test (ANA) Performance Rank") +
  ylim(c(-.3,.17))
child_ameplot
xhist <- axis_canvas(child_ameplot, axis = "x") +
  geom_histogram(data = child_ame, aes(x = mod), binwidth = 5)
child_ameplot <- insert_xaxis_grob(child_ameplot, xhist, position = "bottom")
ggdraw(child_ameplot)


nochild_ameplot <- ggplot(nochild_ame, aes(x = mod, y = ame)) +
  geom_line(size = 2) +
  geom_hline(yintercept = 0, lty = 2, size = 1) +
  geom_ribbon(aes(ymin = ame - 1.96 * ame_se,
                  ymax = ame + 1.96 * ame_se), alpha = .1) +
  geom_pointrange(data = bytercile_nochild, aes(x = median, y = estimate,
                                        ymax = estimate + 1.96 * std.error,
                                        ymin = estimate - 1.96 * std.error),
                  size = .5) +
  theme_minimal() +
  ylab("Average Marginal Effect of Information") +
  xlab("Standardized Test (ANA) Performance Rank") +
  ylim(c(-.3,.17))
xhist <- axis_canvas(nochild_ameplot, axis = "x") +
  geom_histogram(data = nochild_ame, aes(x = mod), binwidth = 5)
nochild_ameplot <- insert_xaxis_grob(nochild_ameplot, xhist, position = "bottom")

#PDF Version
pdf("./figures/ame_child_vs_nochild_plot.pdf", width = 12, height = 5)
plot_grid(child_ameplot, nochild_ameplot,
          labels = c("Child in Local School", "No Child in Local School"),
          align = "h")
dev.off()

#EPS Version
cairo_ps(filename = "./figures/ame_child_vs_nochild_plot.eps",
         width = 12, height = 5, pointsize = 12,
         fallback_resolution = 300)
plot_grid(child_ameplot, nochild_ameplot,
          labels = c("Child in Local School", "No Child in Local School"),
          align = "h")
dev.off()


##Check triple interaction

triple_estimdata <- make_estimdata(ana_contr_all,
                                      dv = "vote", covars = c(covars, "codesetor"))
triple_covars <- lasso_covars(data = triple_estimdata,
                                 dv = "vote",  covars = c(covars, "codesetor"), treat = "ana")
triple_estimdata$rank_child <- triple_estimdata$ana_rank * triple_estimdata$child_school
triple <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + rank_child + child_school + ",
                                                              paste(edu_child_covars, collapse = " + "))),
                    data = triple_estimdata, weights = ipw)

triple_nocovars <- lm_robust(vote ~ ana * demean(ana_rank) * child_school, data = ana_contr_all,
                             fixed_effects = codesetor, weights = ipw)


##Create table
texreg(list(edu, edu_child, edu_nochild, triple),
       include.ci = FALSE,
       custom.coef.map = list("ana" = "Treatment", "ana:ana_rank_c" = "Treatment x Rank",
                              "ana:rank_child_c" = "Treatment x Rank x Parents"),
       include.rmse = FALSE, include.rsquared = FALSE, include.adjrs = FALSE,
       digits = 4, table = FALSE, booktabs = TRUE, use.packages = FALSE, fontsize = "small",
       custom.model.names = c("All", "Parents", "Not Parents", "All"),
       stars = c(.1, .05, .01),
       file = "./tables/exp_mainresults.tex")


## Models without covariate adjustment (for Appendix)

texreg(list(edu_nocovars, edu_child_nocovars, edu_nochild_nocovars, triple_nocovars),
       include.ci = FALSE,
       custom.coef.map = list("ana" = "Treatment", "ana:demean(ana_rank)" = "Treatment x Rank",
                              "ana:demean(ana_rank):child_school" = "Treatment x Rank x Parents"),
       omit.coef  = "(^demean\\(ana_rank\\))|(^child_school)|(^demean\\(ana_rank\\):child_school)",
       include.rmse = FALSE, include.rsquared = FALSE, include.adjrs = FALSE,
       digits = 4, table = FALSE, booktabs = TRUE, use.packages = FALSE, fontsize = "small",
       custom.model.names = c("All", "Parents", "Not Parents", "All"),
       stars = c(.1, .05, .01),
       file = "./tables/exp_mainresults_nocovars.tex"
       )

## Including Other Interactions

edu_otherinter <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + income + years_edu + age + ",
                                                                   paste(edu_covars, collapse = " + "))),
                         data = edu_estimdata, weights = ipw)

edu_child_otherinter <- lm_lin(vote ~ ana,
                               covariates = as.formula(paste("~ ana_rank + income + years_edu + age + ",
                                                             paste(edu_child_covars, collapse = " + "))),
                               data = edu_child_estimdata, weights = ipw)


edu_nochild_otherinter <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + income + years_edu + age + ",
                                                                           paste(edu_nochild_covars, collapse = " + "))),
                                 data = edu_nochild_estimdata, weights = ipw)


texreg(list(edu_otherinter, edu_child_otherinter, edu_nochild_otherinter),
       include.ci = FALSE,
       custom.coef.map = list("ana" = "Treatment", "ana:ana_rank_c" = "Treatment x Rank"),
       include.rmse = FALSE, include.rsquared = FALSE, include.adjrs = FALSE,
       digits = 4, table = FALSE, booktabs = TRUE, use.packages = FALSE, fontsize = "small",
       custom.model.names = c("All", "Parents", "Not Parents"),
       stars = c(.1, .05, .01),
       file = "./tables/exp_mainresults_otherinter.tex"
       )

## For Appendix#####


# Effect by Gap Between Prior and Performance -----------------------------

edu_prior_estimdata <- make_estimdata(ana_contr_all,
                                      dv = "vote", covars = c(covars, "codesetor", "ana_prior_gap"))
edu_prior_estimdata$priorgap_tercile <- ana_contr_all$priorgap_tercile
edu_prior_covars <- lasso_covars(data = edu_prior_estimdata,
                                 dv = "vote",  covars = c(covars, "codesetor", "ana_prior_gap"), treat = "ana")
edu_prior <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_prior_gap +",
                                                              paste(edu_prior_covars, collapse = " + "))),
                    data = edu_prior_estimdata, weights = ipw)
##Divide Sample into Tercile by Ana Prior Gap and estimate effect in each bin
bytercile_gap <- group_by(edu_prior_estimdata, priorgap_tercile) %>%
  nest()
bytercile_gap$median <- map_dbl(bytercile_gap$data, ~ median(.$ana_prior_gap))
bytercile_gap <- bytercile_gap %>%
  mutate(results = map(data, est_func))
bytercile_gap <- unnest(bytercile_gap, results)

get_ame_prior <- function(data, mod){
  ame <- coef(mod)["ana"] + (data$ana_prior_gap - mod$scaled_center["ana_prior_gap"]) * coef(mod)["ana:ana_prior_gap_c"]
  ame_var <- vcov(mod)["ana", "ana"] +
    (data$ana_prior_gap - mod$scaled_center["ana_prior_gap"])^2 * vcov(mod)["ana:ana_prior_gap_c","ana:ana_prior_gap_c"] +
    2 * (data$ana_prior_gap - mod$scaled_center["ana_prior_gap"]) * vcov(mod)["ana:ana_prior_gap_c", "ana"]
  ame_plotdata <- tibble(mod = data$ana_prior_gap, ame = unlist(ame), ame_se = sqrt(unlist(ame_var)))
  ame_plotdata
}

eduprior_ame <- get_ame_prior(edu_prior_estimdata, edu_prior)

eduprior_ameplot <- ggplot(eduprior_ame, aes(x = mod, y = ame)) +
  geom_line(size = 2) +
  geom_hline(yintercept = 0, lty = 2, size = 1) +
  geom_ribbon(aes(ymin = ame - 1.96 * ame_se,
                  ymax = ame + 1.96 * ame_se), alpha = .1) +
  geom_pointrange(data = bytercile_gap, aes(x = median, y = estimate,
                                        ymax = estimate + 1.96 * std.error,
                                        ymin = estimate - 1.96 * std.error),
                  size = .5) +
  theme_minimal() +
  ylab("Average Marginal Effect of Information") +
  xlab("Difference between Prior and Performance")
eduprior_ameplot
xhist <- axis_canvas(eduprior_ameplot, axis = "x") +
  geom_histogram(data = eduprior_ame, aes(x = mod), binwidth = 5)
eduprior_ameplot <- insert_xaxis_grob(eduprior_ameplot, xhist, position = "bottom")
ggdraw(eduprior_ameplot)


pdf("./figures/ame_edugap.pdf", width = 6, height = 4)
ggdraw(eduprior_ameplot)
dev.off()


# Balance -----------------------------------------------------------------


#Remove factor covariates
bal_covars <- covars[covars %in% c("vote_2014", "muni_biggest_prob", "race", "religion") == FALSE]
bal_data <- filter(panel_noattrit,  accounts != 1)
bal_data$codesetor <- as.factor(bal_data$codesetor)


perm_once <- function(covar){
  bal_data <- group_by(bal_data, codesetor) %>%
    mutate(perm_ana = sample(ana))
  bal_data$covar <- bal_data[[covar]]
  lfe::felm(covar ~ perm_ana | codesetor, data = bal_data, weights = bal_data$ipw)$STATS$covar$tval
}

perm_test <- function(covar){
  bal_data$covar <- bal_data[[covar]]
  mod <- lfe::felm(covar ~ ana | codesetor, data = bal_data, weights = bal_data$ipw)
  test_stat <- mod$STATS$covar$tval
  perm_dist <- map_dbl(1:num_perms, ~ perm_once(covar))
  perm_pvalue <- mean(abs(perm_dist) >= abs(test_stat))
  ate_est <-  coef(mod)[["ana"]]
  sd <- sd(bal_data$covar, na.rm = TRUE)
  ate_se <- mod$STATS$covar$se[["ana"]]
  tibble(ate_est = ate_est, ate_se = ate_se, sd = sd, perm_pvalue = perm_pvalue)
}

bal_results <-  future_map_dfr(bal_covars, perm_test, .progress = TRUE)
bal_results$var <- bal_covars


bal_table <- select(bal_results, var, ate_est, sd, ate_se, perm_pvalue) %>%
  arrange(-1 * abs(ate_est/ sd)) %>%
  huxtable()
bal_table <- select(bal_table, var, ate_est, sd, ate_se, perm_pvalue)
align(bal_table[, -1]) <- "center"
bal_table$var <- paste0("\\texttt{", bal_table$var, "}")
bal_table$var <- gsub("_", "\\\\_", bal_table$var)
escape_contents(bal_table[, 1]) <- FALSE
number_format(bal_table[, 1]) <- 0
names(bal_table) <- c("Variable", "Mean Difference", "SD", "SE", "Permutation p-value")
bal_table <- add_colnames(bal_table)
bold(bal_table[1,]) <- TRUE
bottom_border(bal_table[1,]) <- 1
#width(bal_table) <- 1
bal_table <-  bal_table[-2,]
font_size(bal_table) <- 8
bal_table <- set_all_padding(bal_table, 0)

cat(to_latex(bal_table, tabular_only = TRUE), file = "./tables/exp_balance.tex")

#Results for Pre-Registered Hypotheses #####

permute_ana <- function(data){
  group_by(data)%>%
    mutate(ana = sample(ana))
}
perm_mod <- function(mod, data, selected_covars){
  selected_covars <- selected_covars
  future_map_dfr(1:num_perms, ~ tidy(update(mod, data = permute_ana(data))),
                 .id = "permutation", .options = future_options(packages = "estimatr"))
}
perm_mod <- quietly(perm_mod)


##Recipe package function to create estimation dataset
##Remove "ana_tercile" since not pre-specified
make_estimdata <- function(data, dv, covars){
  data <- data[, names(data) %in% c(covars, dv, "ana", "ipw", "ana_rank")]
  rec <- recipe(as.formula(paste(dv, " ~ ", paste(c(covars), collapse = "+"))),
                data = data) %>%
    step_dummy(all_nominal()) %>%
    step_meanimpute(all_predictors()) %>%
    step_zv(all_predictors()) %>%
    prep()
  estim_data <- juice(rec)
  estim_data$ana <- data$ana
  estim_data$ana_rank <- data$ana_rank
  estim_data$ipw <- data$ipw
  estim_data
}

num_perms <- 10000


#Family 1
##Hypothesis 1a
hyp1a_data <- make_estimdata(ana_contr_good, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp1a_data, dv = "vote",
                             covars = c(covars, "codesetor"), treat = "ana")

hyp1a <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank +",
                                                        paste(selected_covars, collapse = " + "))),
                data = hyp1a_data, weights = ipw)
hyp1a_permdist <- perm_mod(hyp1a, hyp1a_data, selected_covars)$result


##Hypothesis 1b
hyp1b_data <- make_estimdata(ana_contr_bad, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp1b_data, dv = "vote",
                             covars = c(covars, "codesetor"), treat = "ana")

hyp1b <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank +",
                                                          paste(selected_covars, collapse = " + "))),
                data = hyp1b_data, weights = ipw)
hyp1b_permdist <- perm_mod(hyp1b, hyp1b_data, selected_covars)$result

#Family 2
##Hypothesis 2
###Good News
hyp2good_data <- make_estimdata(ana_contr_good, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp2good_data, dv = "vote",
                             covars = c(covars, "codesetor"), treat = "ana")

hyp2good <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + uncertain_edu_baseline + ",
                                                          paste(selected_covars, collapse = " + "))),
                data = hyp2good_data, weights = ipw)
hyp2good_permdist <- perm_mod(hyp2good, hyp2good_data, selected_covars)$result

###Bad News
hyp2bad_data <- make_estimdata(ana_contr_bad, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp2bad_data, dv = "vote",
                                covars = c(covars, "codesetor"), treat = "ana")
hyp2bad <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + uncertain_edu_baseline + ",
                                                             paste(selected_covars, collapse = " + "))),
                   data = hyp2bad_data, weights = ipw)
hyp2bad_permdist <- perm_mod(hyp2bad, hyp2bad_data, selected_covars)$result


##Hypothesis 3
###Good News
hyp3good_data <- make_estimdata(ana_contr_good, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp3good_data, dv = "vote",
                             covars = c(covars, "codesetor"), treat = "ana")
hyp3good <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + edu_rank1 + ",
                                                          paste(selected_covars, collapse = " + "))),
                data = hyp3good_data, weights = ipw)
hyp3good_permdist <- perm_mod(hyp3good, hyp3good_data, selected_covars)$result

###Bad News
hyp3bad_data <- make_estimdata(ana_contr_bad, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp3bad_data, dv = "vote",
                                covars = c(covars, "codesetor"), treat = "ana")
hyp3bad <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + edu_rank1 + ",
                                                             paste(selected_covars, collapse = " + "))),
                   data = hyp3bad_data, weights = ipw)
hyp3bad_permdist <- perm_mod(hyp3bad, hyp3bad_data, selected_covars)$result

##Hypothesis 4
###Good News
hyp4good_data <- make_estimdata(ana_contr_good, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp4good_data, dv = "vote",
                                covars = c(covars, "codesetor"), treat = "ana")
hyp4good_data$edu_rank1_uncertain_edu <- with(hyp4good_data, edu_rank1 * uncertain_edu_baseline)
hyp4good <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + edu_rank1 + uncertain_edu_baseline + edu_rank1_uncertain_edu + ",
                                                             paste(selected_covars, collapse = " + "))),
                data = hyp4good_data, weights = ipw)
hyp4good_permdist <- perm_mod(hyp4good, hyp4good_data, selected_covars)$result

###Bad News
hyp4bad_data <- make_estimdata(ana_contr_bad, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp4bad_data, dv = "vote",
                                covars = c(covars, "codesetor"), treat = "ana")
hyp4bad_data$edu_rank1_uncertain_edu <- with(hyp4bad_data, edu_rank1 * uncertain_edu_baseline)
hyp4bad <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + edu_rank1 + uncertain_edu_baseline + edu_rank1_uncertain_edu + ",
                                                             paste(selected_covars, collapse = " + "))),
                   data = hyp4bad_data, weights = ipw)
hyp4bad_permdist <- perm_mod(hyp4bad, hyp4bad_data, selected_covars)$result

##Hypothesis 5a
###Good News
hyp5agood_data <- make_estimdata(ana_contr_good, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp5agood_data, dv = "vote",
                                covars = c(covars, "codesetor"), treat = "ana")
hyp5agood <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + child_school + ",
                                                             paste(selected_covars, collapse = " + "))),
                   data = hyp5agood_data, weights = ipw)
hyp5agood_permdist <- perm_mod(hyp5agood, hyp5agood_data, selected_covars)$result

###Bad News
hyp5abad_data <- make_estimdata(ana_contr_bad, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp5abad_data, dv = "vote",
                                covars = c(covars, "codesetor"), treat = "ana")
hyp5abad <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + child_school + ",
                                                             paste(selected_covars, collapse = " + "))),
                   data = hyp5abad_data, weights = ipw)
hyp5abad_permdist <- perm_mod(hyp5abad, hyp5abad_data, selected_covars)$result

##Hypothesis 5b
###Good News
hyp5bgood_data <- make_estimdata(ana_contr_good, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp5bgood_data, dv = "vote",
                                 covars = c(covars, "codesetor"), treat = "ana")
hyp5bgood_data$child_school_uncertain_edu <- with(hyp5bgood_data, child_school * uncertain_edu_baseline)
hyp5bgood <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + child_school + uncertain_edu_baseline + child_school_uncertain_edu + ",
                                                              paste(selected_covars, collapse = " + "))),
                    data = hyp5bgood_data, weights = ipw)
hyp5bgood_permdist <- perm_mod(hyp5bgood, hyp5bgood_data, selected_covars)$result

###Bad News
hyp5bbad_data <- make_estimdata(ana_contr_bad, dv = "vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp5bbad_data, dv = "vote",
                                 covars = c(covars, "codesetor"), treat = "ana")
hyp5bbad_data$child_school_uncertain_edu <- with(hyp5bbad_data, child_school * uncertain_edu_baseline)
hyp5bbad <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + child_school + uncertain_edu_baseline + child_school_uncertain_edu + ",
                                                              paste(selected_covars, collapse = " + "))),
                    data = hyp5bbad_data, weights = ipw)
hyp5bbad_permdist <- perm_mod(hyp5bbad, hyp5bbad_data, selected_covars)$result

#Family 3
##Hypothesis 7a
hyp7a_data <- make_estimdata(ana_contr_good, dv = "turnout", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp7a_data, dv = "turnout",
                             covars = c(covars, "codesetor"), treat = "ana")
hyp7a <- lm_lin(turnout ~ ana, covariates = as.formula(paste("~ ana_rank +",
                                                          paste(selected_covars, collapse = " + "))),
                data = hyp7a_data, weights = ipw, se_type = "HC1")
hyp7a_permdist <- perm_mod(hyp7a, hyp7a_data, selected_covars)$result

##Hypothesis 7b
hyp7b_data <- make_estimdata(ana_contr_bad, dv = "turnout", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp7b_data, dv = "turnout",
                             covars = c(covars, "codesetor"), treat = "ana")
hyp7b <- lm_lin(turnout ~ ana, covariates = as.formula(paste("~ ana_rank +",
                                                          paste(selected_covars, collapse = " + "))),
                data = hyp7b_data, weights = ipw, se_type = "HC1")
##Use HC1 because HC2 SE generates NA
hyp7b_permdist <- perm_mod(hyp7b, hyp7b_data, selected_covars)$result

#Family 4
##Hypothesis 8a
hyp8a_data <- make_estimdata(ana_contr_bad, dv = "valid_vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp8a_data, dv = "valid_vote",
                             covars = c(covars, "codesetor"), treat = "ana")
hyp8a <- lm_lin(valid_vote ~ ana, covariates = as.formula(paste("~ ana_rank +",
                                                             paste(selected_covars, collapse = " + "))),
                data = hyp8a_data, weights = ipw, se_type = "HC1")
hyp8a_permdist <- perm_mod(hyp8a, hyp8a_data, selected_covars)$result

##Hypothesis 8b
hyp8b_data <- make_estimdata(ana_contr_bad, dv = "valid_vote", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp8b_data, dv = "valid_vote",
                             covars = c(covars, "codesetor"), treat = "ana")
hyp8b <- lm_lin(valid_vote ~ ana, covariates = as.formula(paste("~ ana_rank + govt_eval_baseline + ",
                                                                paste(selected_covars, collapse = " + "))),
                data = hyp8b_data, weights = ipw, se_type = "HC1")
hyp8b_permdist <- perm_mod(hyp8b, hyp8b_data, selected_covars)$result

##Hypothesis 9
###Good News
hyp9good_data <- make_estimdata(ana_contr_good, dv = "edu_eval", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp9good_data, dv = "edu_eval",
                             covars = c(covars, "codesetor"), treat = "ana")

hyp9good <- lm_lin(edu_eval ~ ana, covariates = as.formula(paste("~ ana_rank + ",
                                                          paste(selected_covars, collapse = " + "))),
                data = hyp9good_data, weights = ipw, se_type = "HC1")
hyp9good_permdist <- perm_mod(hyp9good, hyp9good_data, selected_covars)$result

###Bad News
hyp9bad_data <- make_estimdata(ana_contr_bad, dv = "edu_eval", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp9bad_data, dv = "edu_eval",
                                covars = c(covars, "codesetor"), treat = "ana")
hyp9bad <- lm_lin(edu_eval ~ ana, covariates = as.formula(paste("~ ana_rank + ",
                                                             paste(selected_covars, collapse = " + "))),
                   data = hyp9bad_data, weights = ipw)
hyp9bad_permdist <- perm_mod(hyp9bad, hyp9bad_data, selected_covars)$result

##Hypothesis 10
###Good News
hyp10good_data <- make_estimdata(ana_contr_good, dv = "edu_eval", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp10good_data, dv = "edu_eval",
                             covars = c(covars, "codesetor"), treat = "ana")

hyp10good <- lm_lin(edu_eval ~ ana, covariates = as.formula(paste("~ ana_rank + uncertain_edu_baseline + ",
                                                          paste(selected_covars, collapse = " + "))),
                data = hyp10good_data, weights = ipw, se_type = "HC1")
hyp10good_permdist <- perm_mod(hyp10good, hyp10good_data, selected_covars)$result

###Bad News
hyp10bad_data <- make_estimdata(ana_contr_bad, dv = "edu_eval", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp10bad_data, dv = "edu_eval",
                                covars = c(covars, "codesetor"), treat = "ana")
hyp10bad <- lm_lin(edu_eval ~ ana, covariates = as.formula(paste("~ ana_rank + uncertain_edu_baseline + ",
                                                             paste(selected_covars, collapse = " + "))),
                   data = hyp10bad_data, weights = ipw)
hyp10bad_permdist <- perm_mod(hyp10bad, hyp10bad_data, selected_covars)$result

##Hypothesis 11
hyp11_data <- make_estimdata(ana_contr_all, dv = "uncertain_edu_eval", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp11_data, dv = "uncertain_edu_eval",
                                covars = c(covars, "codesetor"), treat = "ana")
hyp11 <- lm_lin(uncertain_edu_eval ~ ana, covariates = as.formula(paste("~ ana_rank + ",
                                                                 paste(selected_covars, collapse = " + "))),
                   data = hyp11_data, weights = ipw, se_type = "HC1")
hyp11_permdist <- perm_mod(hyp11, hyp11_data, selected_covars)$result

##Hypothesis 12
hyp12_data <- make_estimdata(ana_contr_all, dv = "info_import", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp12_data, dv = "info_import",
                             covars = c(covars, "codesetor"), treat = "ana")
hyp12 <- lm_lin(info_import ~ ana, covariates = as.formula(paste("~ ana_rank + ",
                                                                        paste(selected_covars, collapse = " + "))),
                data = hyp12_data, weights = ipw, se_type = "HC1")
hyp12_permdist <- perm_mod(hyp12, hyp12_data, selected_covars)$result

##Hypothesis 12
hyp13_data <- make_estimdata(ana_contr_all, dv = "ana_correct", covars = c(covars, "codesetor"))
selected_covars <- lasso_covars(data = hyp13_data, dv = "ana_correct",
                             covars = c(covars, "codesetor"), treat = "ana")
hyp13 <- lm_lin(ana_correct ~ ana, covariates = as.formula(paste("~ ana_rank + ",
                                                                 paste(selected_covars, collapse = " + "))),
                data = hyp13_data, weights = ipw, se_type = "HC1")
hyp13_permdist <- perm_mod(hyp13, hyp13_data, selected_covars)$result


#Create table

models <- list(hyp1a, hyp1b, hyp2good, hyp2bad, hyp3good, hyp3bad,
               hyp4good, hyp4bad, hyp5agood, hyp5abad, hyp5bgood, hyp5bbad,
               hyp7a, hyp7b, hyp8a, hyp8b, hyp9good, hyp9bad,
               hyp10good, hyp10bad, hyp11, hyp12, hyp13)
perm_dists_allvars <- list(hyp1a_permdist, hyp1b_permdist, hyp2good_permdist,
                           hyp2bad_permdist, hyp3good_permdist, hyp3bad_permdist,
                           hyp4good_permdist, hyp4bad_permdist, hyp5agood_permdist,
                           hyp5abad_permdist, hyp5bgood_permdist, hyp5bbad_permdist,
                           hyp7a_permdist, hyp7b_permdist, hyp8a_permdist, hyp8b_permdist,
                           hyp9good_permdist, hyp9bad_permdist,
                           hyp10good_permdist, hyp10bad_permdist, hyp11_permdist,
                           hyp12_permdist, hyp13_permdist)
table_vars <- c("ana", "ana", "ana:ana_rank_c", "ana:ana_rank_c",
                "ana:edu_rank1_c", "ana:edu_rank1_c", "ana:edu_rank1_uncertain_edu_c",
                "ana:edu_rank1_uncertain_edu_c", "ana:child_school_c", "ana:child_school_c",
                "ana:child_school_uncertain_edu_c", "ana:child_school_uncertain_edu_c",
                "ana", "ana", "ana", "ana:govt_eval_baseline_c", "ana", "ana",
                "ana:uncertain_edu_baseline_c", "ana:uncertain_edu_baseline_c",
                "ana", "ana", "ana")
table_ests <- map2_dbl(models, table_vars, ~ .x$coefficients[.y])
table_se <- map2_dbl(models, table_vars, ~ .x$std.error[.y])
table_t <- map2_dbl(models, table_vars, ~ .x$statistic[.y])
table_p <- map2_dbl(models, table_vars, ~ .x$p.value[.y])

perm_dists <- map2(perm_dists_allvars,
                   table_vars, ~ filter(.x, term == .y)$statistic)
table_perm_p <- map2_dbl(table_t, perm_dists, ~ sum(abs(.y) >= abs(.x))/(length(.y) + 1))


prereg_table <- tibble(Family = c(1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4),
       Hypothesis = c("1a", "1b", "2", "2", "3", "3", "4", "4", "5a", "5a", "5b", "5b",
                      "7a", "7b", "8a", "8b", "9", "9", "10", "10", "11", "12", "13"),
       Outcome = c("vote", "vote", "vote", "vote", "vote", "vote", "vote", "vote",
                   "vote", "vote", "vote", "vote", "turnout", "turnout", "valid_vote", "valid_vote",
                   "edu_eval", "edu_eval", "edu_eval", "edu_eval", "uncertain_edu_eval", "info_import",
                   "ana_correct"),
       Parameter = c("ana", "ana", "ana x uncertain_edu", "ana x uncertain_edu",
                     "ana x edu_rank1", "ana x edu_rank1", "ana x edu_rank1 x uncertain_edu",
                     "ana x edu_rank1 x uncertain_edu", "ana x child_school", "ana x child_school",
                     "ana x child_school x uncertain_edu", "ana x child_school x uncertain_edu",
                     "ana", "ana", "ana", "ana x govt_eval", "ana", "ana",
                     "ana x uncertain_edu", "ana x uncertain_edu", "ana", "ana", "ana"),
       Subset = c("Good News", "Bad News", "Good News", "Bad News", "Good News",
                  "Bad News", "Good News", "Bad News", "Good News", "Bad News",
                  "Good News", "Bad News", "Good News", "Bad News", "Bad News", "Bad News",
                  "Good News", "Bad News", "Good News", "Bad News", "All", "All", "All"),
       Estimate = table_ests,
       SE = table_se,
       `p-value` = table_p) %>%
  mutate(Outcome = cell_spec(Outcome, "latex", monospace = TRUE),
         Parameter = cell_spec(Parameter, "latex", monospace = TRUE)) %>%
  kable("latex", booktabs = TRUE,
        digits = 3, escape = FALSE,
        linesep = "")


cat(prereg_table, file = "./tables/prereg_results.tex")


# Attrition Analysis ------------------------------------------------------

edu_attrit <- filter(panel_noattrit, accounts != 1)
#Impute missing covariate data
covar_data_attrit <- mlr::impute(edu_attrit[,covars],
                          classes = list(integer = mlr::imputeMean(), numeric = mlr::imputeMean(),
                                         factor = mlr::imputeMode(), character = mlr::imputeMode()))$data[, covars]
edu_attrit <- bind_cols(edu_attrit[, ! names(edu_attrit) %in% covars], covar_data_attrit)

attrit <- lm_robust(attrited ~ ana,
                          data = edu_attrit, weights = ipw,
                          fixed_effects = codesetor)
attrit_inter <- lm_robust(attrited ~ ana * demean(ana_rank) + ana * demean(vote_2012) +
                            ana * demean(govt_eval_baseline) + ana * demean(child_school) +
                            ana * demean(confid_muni),
                          data = edu_attrit, weights = ipw,
                          fixed_effects = codesetor)

attrit_table <- huxreg(attrit, attrit_inter,
       coefs = c("Treatment" = "ana",
                 "Treatment x rank" = "ana:demean(ana_rank)",
                 "Treatment x vote_2012" = "ana:demean(vote_2012)",
                 "Treatment x govt_eval_baseline" = "ana:demean(govt_eval_baseline)",
                 "Treatment x child_school" = "ana:demean(child_school)",
                 "Treatment x confid_muni" = "ana:demean(confid_muni)"),
       statistics = c("N" = "nobs"),
       stars = c(`***` = 0.01, `**` = 0.05, `*` = 0.1))

cat(to_latex(attrit_table, tabular_only = TRUE), file = "./tables/exp_attrition.tex")


## Ceiling Effects Analysis ------------------------------------------------------
load("./data/ana.RData")

ana$ana_baseline <- ana$math13.avg+ana$reading13.avg
ana_baseline <- select(ana, muni, ana_baseline) %>%
  arrange(muni) %>%
  mutate(muni = as.character(muni))
ana_baseline$muni[ana_baseline$muni == "LAGOA DO GATOS"] <- "LAGOA DOS GATOS"

exp_ana_baseline <- select(panel, muni, codeibge) %>%
  mutate(muni = stringi::stri_trans_general(toupper(muni), "latin-ascii")) %>%
  distinct() %>%
  left_join(ana_baseline) %>%
  select(-muni)

ana_contr_all <- left_join(ana_contr_all, exp_ana_baseline)


edu_anabaseline_estimdata <- make_estimdata(ana_contr_all, dv = "vote", covars = c(covars, "codesetor", "ana_baseline"))
edu_anabaseline_covars <- lasso_covars(data = edu_anabaseline_estimdata, dv = "vote",
                           covars = c(covars, "codesetor", "ana_baseline"), treat = "ana")
edu_anabaseline <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + ana_baseline +",
                                                        paste(edu_covars, collapse = " + "))),
              data = edu_anabaseline_estimdata, weights = ipw)


edu_child_anabaseline_estimdata <- make_estimdata(filter(ana_contr_all, child_school ==1),
                                      dv = "vote", covars = c(covars, "codesetor", "ana_baseline"))
edu_child_anabaseline_covars <- lasso_covars(data = edu_child_anabaseline_estimdata,
                                 dv = "vote",  covars = c(covars, "codesetor", "ana_baseline"), treat = "ana")
edu_child_anabaseline <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + ana_baseline + ",
                                                        paste(edu_child_anabaseline_covars, collapse = " + "))),
              data = edu_child_anabaseline_estimdata, weights = ipw)


edu_nochild_anabaseline_estimdata <- make_estimdata(filter(ana_contr_all, child_school == 0),
                                      dv = "vote", covars = c(covars, "codesetor", "ana_baseline"))
edu_nochild_anabaseline_covars <- lasso_covars(data = edu_nochild_anabaseline_estimdata,
                                 dv = "vote",  covars = c(covars, "codesetor", "anabaseline"), treat = "ana")
edu_nochild_anabaseline <- lm_lin(vote ~ ana, covariates = as.formula(paste("~ ana_rank + ana_baseline + ",
                                                              paste(edu_nochild_anabaseline_covars, collapse = " + "))),
                    data = edu_nochild_anabaseline_estimdata, weights = ipw)


texreg(list(edu_anabaseline, edu_child_anabaseline, edu_nochild_anabaseline),
       include.ci = FALSE,
       custom.coef.map = list("ana" = "Treatment", "ana:ana_rank_c" = "Treatment x Rank"),
       include.rmse = FALSE, include.rsquared = FALSE, include.adjrs = FALSE,
       digits = 4, table = FALSE, booktabs = TRUE, use.packages = FALSE, fontsize = "small",
       custom.model.names = c("All", "Parents", "Not Parents"),
       stars = c(.1, .05, .01),
       file = "./tables/exp_mainresults_anabaseline.tex")


# R version 3.6.3 (2020-02-29)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Ubuntu 18.04.4 LTS
#
# Matrix products: default
# BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
# LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
#
# locale:
#   [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
# [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C
# [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#
# attached base packages:
#   [1] stats     graphics  grDevices utils     datasets  methods   base
#
# other attached packages:
#   [1] furrr_0.1.0        future_1.17.0      yardstick_0.0.6    rsample_0.0.7      recipes_0.1.13
# [6] parsnip_0.1.1.9000 kableExtra_1.1.0   texreg_1.37.5      estimatr_0.22.0    huxtable_5.0.0
# [11] cowplot_1.0.0      forcats_0.5.0      stringr_1.4.0      dplyr_1.0.0        purrr_0.3.4
# [16] readr_1.3.1        tidyr_1.1.0        tibble_3.0.1       ggplot2_3.3.2      tidyverse_1.3.0
#
# loaded via a namespace (and not attached):
#   [1] colorspace_1.4-1   ellipsis_0.3.1     class_7.3-17       fs_1.4.1           rstudioapi_0.11
# [6] listenv_0.8.0      farver_2.0.3       ParamHelpers_1.14  prodlim_2019.11.13 fansi_0.4.1
# [11] lubridate_1.7.9    xml2_1.3.2         codetools_0.2-16   splines_3.6.3      knitr_1.29
# [16] Formula_1.2-3      jsonlite_1.6.1     pROC_1.16.2        packrat_0.5.0      broom_0.5.6
# [21] dbplyr_1.4.4       compiler_3.6.3     httr_1.4.1         backports_1.1.8    assertthat_0.2.1
# [26] Matrix_1.2-18      lfe_2.8-5          cli_2.0.2          htmltools_0.5.0    tools_3.6.3
# [31] gtable_0.3.0       glue_1.4.1         fastmatch_1.1-0    Rcpp_1.0.4.6       parallelMap_1.5.0
# [36] cellranger_1.1.0   vctrs_0.3.1        nlme_3.1-147       iterators_1.0.12   timeDate_3043.102
# [41] gower_0.2.1        xfun_0.15          mlr_2.17.1         globals_0.12.5     rvest_0.3.5
# [46] lifecycle_0.2.0    MASS_7.3-51.6      zoo_1.8-8          scales_1.1.1       ipred_0.9-9
# [51] hms_0.5.3          parallel_3.6.3     sandwich_2.5-1     BBmisc_1.11        rpart_4.1-15
# [56] stringi_1.4.6      foreach_1.5.0      checkmate_2.0.0    lava_1.6.7         shape_1.4.4
# [61] rlang_0.4.6        pkgconfig_2.0.3    commonmark_1.7     evaluate_0.14      lattice_0.20-41
# [66] labeling_0.3       tidyselect_1.1.0   plyr_1.8.6         magrittr_1.5       R6_2.4.1
# [71] generics_0.0.2     DBI_1.1.0          pillar_1.4.4       haven_2.3.1        withr_2.2.0
# [76] survival_3.2-3     nnet_7.3-14        modelr_0.1.8       crayon_1.3.4       rmarkdown_2.3
# [81] grid_3.6.3         readxl_1.3.1       data.table_1.12.8  blob_1.2.1         reprex_0.3.0
# [86] digest_0.6.25      webshot_0.5.2      xtable_1.8-4       munsell_0.5.0      glmnet_4.0-2
# [91] viridisLite_0.3.0
# >
